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Abstract 

High-order  non-reflecting  boundary  conditions  are  introduced  to  create  a  finite  computational  space  and 
for  the  solution  of  dispersive  waves  using  a  spectral  element  formulation  with  high-order  time  integration. 
Numerical  examples  are  used  to  demonstrate  the  synergy  of  using  high-order  spatial,  time,  and  boundary 
discretization.  We  show  that  by  balancing  all  numerical  errors  involved,  high-order  accuracy  can  be  achieved 
for  unbounded  domain  problems  in  polar  coordinate  systems. 
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1.  Introduction 

The  numerical  solution  of  a  wave  propagation  problem  in  a  very  large  or  unbounded  domain  provides  a 
challenging  computational  difficulty  -  namely,  solving  the  problem  on  a  finite  computational  domain  while 
maintaining  the  true  essence  of  the  solution.  One  of  the  modern  techniques  that  has  garnered  a  significant 
amount  of  attention  in  handling  this  challenge  is  the  absorbing  or  non-reflecting  boundary  condition  (NRBC) 
method.  In  using  this  method,  the  original  infinite  domain  is  truncated  by  an  artificial  boundary  B ,  resulting 
in  a  finite  computational  domain  Q  and  the  residual  domain  D. 

When  truncating  the  domain,  the  modeler  must  devise  boundary  conditions  for  the  truncated  domain. 
Of  course,  by  imposing  a  boundary  where  one  does  not  physically  exist,  the  problem  is  changed  -  and  unless 
chosen  carefully,  would  certainly  be  expected  to  pollute  the  solution  as  the  problem  evolves  and  impinges 
on  the  boundary.  For  this  reason,  much  effort  has  and  continues  to  be  exerted  on  finding  stable,  efficient, 
accurate  and  practical  means  of  reducing  this  reflection  through  so-called  NRBCs  [1]. 

Several  high- order  NRBCs  have  been  devised  to  reduce  spurious  reflections  that  would  pollute  the  so¬ 
lution.  Beginning  in  the  late  1980’s,  the  well-known  Engquist-Majda  [2]  and  Bayliss-Turkel  conditions  [3] 
gave  way  to  Collino’s  [4]  low  derivative,  auxiliary  variable  formulation  for  the  2D  scalar  wave  equation.  This 
sparked  a  flurry  of  activity  in  an  effort  to  find  quality,  high-order  NRBCs  that  were  easily  implementable. 
See  [5]- [6]  for  reviews  on  the  subject.  See  also  Higdon’s  papers  [7]- [8],  Givoli-Neta  [9],  [10]  and  Neta  et  al 
[11]- 

The  Givoli-Neta  (G-N)  auxiliary  formulation  of  the  Higdon  NRBC  was  implemented  on  a  reduced  form 
of  the  linearized  shallow  water  equations  (SWE)  under  non-zero  advection.  This  formulation  has  been 
previously  demonstrated  in  a  finite  difference  formulation  to  arbitrarily  high  NRBC  order  [12],  however, 
accuracy  gains  realized  by  increasing  the  NRBC  ceased  after  order  2.  The  formulation  used  by  Lindquist  et 
al  [13]  remedied  this  limitation  by  using  a  high-order  treatment  of  space  (SE)  and  time  (Runge-Kutta)  to 
show  the  benefits  of  using  the  lrigh-order  boundary  (G-N)  scheme.  Specifically,  the  interior  and  boundary 
formulations  are  discretized  using  high-order  basis  functions  in  a  stable,  equal-order  interpolation  scheme 
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for  all  the  variables  (this  does  not  violate  the  inf-sup  condition).  High-order  time  integration  is  performed  as 
well  in  an  effort  to  balance  all  of  the  errors  involved  with  the  numerical  solution.  The  computational  effort 
associated  with  the  high-order  boundary  scheme  can  be  shown  to  grow  only  linearly  with  the  order  [14]. 

It  should  be  noted  that  the  only  other  spectral  element,  high-order  boundary  approach  are  Kucherov 
and  Givoli  [15]  and  Lindquist  et  al.  [16]  and  [13].  Kucherov  and  Givoli  demonstrate  exponential  error 
convergence  of  the  classical  wave  equation  on  a  semi-infinite  channel  when  solved  using  spectral  elements 
and  high  order  boundary  treatment  (using  the  H-W  boundary  scheme).  They  further  show  how  the  spectral 
element  formulation  allows  the  NRBC  to  realize  its  true  potential,  prior  masked  by  low  order  numerical 
schemes.  They  note  that,  “Although  it  is  generally  felt  that  there  is  no  need  to  treat  the  time  domain 
‘spectrally’  like  the  spatial  domain,  [they]  feel  that  a  consistently  high-order  treatment  requires  that  the 
entire  approximation  be  spectral,  i.e.  the  convergence  of  all  three  types  of  error  -  the  spatial  and  temporal 
discretization  errors  and  the  ABC  error-  be  exponential.”  [15] 

Subsequent  spectral  element  and  boundary  work  by  the  current  authors  [16]  using  the  G-N  formulation 
of  the  dispersive  wave  equation  on  a  semi-infinite  channel  showed  similar  results  to  those  presented  by 
Kucherov  and  Givoli  and  how  a  high-order  treatment  of  the  time  domain  (up  to  order  10)  produces  additional 
improvements.  These  improvements,  however,  have  their  limits  thus  confirming  the  hypothesis  of  Kucherov 
and  Givoli.  The  key  difference  in  that  work  is  that  we  extend  the  high  order  space,  boundary  and  time 
integration  results  previously  demonstrated  in  a  non-zero  advection  setting  to  one  where  the  wave  medium 
is  not  at  rest. 

In  the  case  of  cylindrical  or  spherical  coordinate  systems,  one  cannot  use  the  auxiliary  functions  suggested 
for  the  cartesian  coordinate  system.  The  reason  is  that  only  in  the  latter  the  coefficients  are  constant.  To 
overcome  this,  van  Joolen  at  al  [17]  extended  the  Hagstrom-Hariharan  (HH)  conditions  to  dispersive  media. 
Here  we  show  how  to  implement  this  idea.  The  following  is  the  outline  of  the  rest  of  this  paper.  In  Section 
2,  the  problem  under  investigation  is  stated.  In  Section  3,  an  overview  of  the  auxiliary  formulation  is 
presented.  In  Section  4,  a  SE  semi-discrete  formulation  that  incorporates  the  NRBC  with  any  desired  order 
is  constructed.  In  Section  5  the  time-integrator  used  to  march  the  equations  in  time  is  discussed.  The 
performance  of  the  method  is  demonstrated  in  Section  6  via  numerical  examples. 


2.  Statement  of  the  Problem 

To  motivate  the  problem  under  consideration,  consider  the  SWE: 

dt.u  +  udxu  +  vdyU  —  fv  =  —g  dxh , 

dtv  +  udxv  +  vdyV  +  fu  =  —g  dyh ,  (1) 

dth  +  dx(Hu)  +  dy(Hv)  =  0. 

We  use  the  following  shorthand  for  partial  derivatives 

dl 

dl  =  — 
a  dai 

The  shallow  water  model  in  its  current  form  is  non-linear.  We  have  three  state  variables:  u(x,  y1 1)  and 

v(x,y,t)  are  the  unknown  velocities  in  the  x  and  y  directions  and  h(x,y,t)  is  the  water  depth  above  a 

reference  value.  Further,  H  is  the  water  depth  as  shown  in  Figure  1  such  that  H  =  hs  +  h,  f  is  the  Coriolis 
parameter,  and  g  is  the  gravity  acceleration. 

Now,  suppose  that  the  bottom  topography  is  flat  such  that  hs  is  constant  and  u  and  v  can  be  described 
by  a  constant  mean  term  and  a  small  0(5)  deviation  from  that  value,  i.e. 

u  =  U  +  u*  v  =  V  +  v*. 

To  be  clear,  U  and  V  are  the  mean  velocities  with  respect  to  the  coordinate  axes.  Using  these  substitutions 
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and  neglecting  any  0(S2)  terms  results  in  the  linearized  form  of  the  SWE: 

dtu*  +  Udxu*  +  VdyU*  -  f(V  +  v*)  =  - g  dxh, 

dtv*  +  Udxv*  +  VdyV*  +  f(U  +  u*)  =  - g  dvh ,  (2) 

dth  +  Udxh  +  Vdyh  +  hs  ( dxu *  +  dyv*)  =  0. 

van  Joolen  shows  in  [18],  how,  using  the  operator: 

D  _  d  d  d 
Dt  dt  ^  dx  ^  dy 

the  SWE  can  be  reduced  to  a  single  variable  Klein  Gordon  equation  (KGE)  equivalent  for  the  wave  per¬ 
turbation  h  under  non-zero,  constant  advection  velocities  U  and  V.  Once  h  has  been  found,  the  unknown 
velocities  can  be  computed  in  a  similar  manner  (see  Pedlosky  [19]  for  details).  The  current  paper  seeks  to 
numerically  solve  this  two-dimensional  advective,  dispersive  wave  equation: 

d2th  +  (U2  -  Cg)  d2h  +  (V2  -  Cg)  d2yh  +  2 Ud2xth  +  2 Vd2yth  +  2UVd2xyh  +  f2h  =  0  (3) 

where  Cg  =  gH ,  using  continuous  Galerkin  methods. 

While  the  assumptions  for  this  reduced  form  of  the  SWE  may  undermine  the  predictive  capability  of  the 
already  simplified  equations  of  fluid  motion,  the  advective  KGE  serves  as  an  important  test-bed  to  extend 
and  improve  non- reflecting  boundary  conditions  for  wave  propagation  in  a  non-static  environment.  When 
shown  to  improve  performance  in  test  cases  as  those  presented  here,  these  conditions  will  then  be  extended 
to  the  full  linearized  SWE  system  and  further  to  include  non-linear  effects. 


3.  Hagstrom  and  Hariharan  auxiliary  variable  formulation 


The  boundary  condition  devised  by  Hagstrom  and  Hariharan  provides  a  systematic  approach  for  con¬ 
structing  boundary  conditions  for  the  standard  two-dimensional  wave  equation.  The  condition  is  based  on 
the  asymptotic  series  representation  (which  does  not  converge  at  any  fixed  radius)  for  an  outgoing  solution 
of  the  wave  equation  (in  polar  coordinates) 


1  d2h  d2h  1  dh  1  d2h 
Cq  dt2  dr 2  r  dr  +  r2  d02 


(4) 


Since  the  boundary  condition  is  asymptotic  by  nature,  valid  for  large  radial  distances  -  this  implies  that 
larger  radial  distances  should  provide  better  NRBC  convergence.  Thompson  et  al.  make  the  observation 
that  “. . .  for  practical  problems,  truncating  the  asymptotic  expansion  after  [J]  terms  provides  solutions 
with  errors  well  below  that  of  the  discretization  error”  [20].  Here,  we  seek  to  significantly  reduce  the 
discretization  error  by  employing  spectral  elements  to  find  the  true  error  convergence  properties  of  the 
NRBC.  In  developing  the  boundary  condition,  Hagstrom  and  Hariharan  construct  a  sequence  of  operators 
that  approximately  annihilate  the  residual  of  the  preceding  element  in  the  sequence,  viewed  as  a  function 
on  the  artificial  boundary.  The  sequence  begins  with  a  first-order  Bayliss-Turkcl  operator  discussed  in  [3]. 
The  boundary  condition  takes  the  form: 


where 


dh 

dr 


=  <l>i 


_  i  d<t>j 
^+1~co  ~dt 


1  dh 
cq  dt 

j 


2  r 


(J^)2 

4r2 


Pj-i 


1 

4r2  dff2 


j  =  1, . . . ,  J  -  1 


4>  o  =  2  h  and  <f>j  =  0. 


(5) 

(6) 


At  first  glance,  this  boundary  formulation  suggests  that  we  should  develop  a  “new”  spectral  element 
formulation  for  the  wave  equation  cast  in  polar  coordinates.  If  we  did  this,  however,  we  would  then  require 
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a  polar  grid  that  would  introduce  additional  complications  such  as  the  method  of  dealing  with  the  degenerate 
quadrilaterals  that  inevitably  occur  at  the  center  of  the  grid.  Of  course  there  are  ways  to  overcome  these 
obstacles,  but  it  would  be  much  more  convenient  to  cast  the  problem  in  the  same  framework  already 
developed.  In  other  words,  we  seek  to  implement  this  boundary  condition  (presented  in  polar  form)  in  our 
unstructured  quadrilateral  formulation  of  the  wave  equation  (in  Cartesian  form). 


4.  Spectral  Element  Method 

The  spectral  element  (SE)  method  is  a  generalized  high-order  finite  element  method  where  the  integration 
and  interpolation  points  are  selected  carefully  in  order  to  yield  accurate  but  efficient  solutions.  For  this 
problem,  we  will  discuss  two  formulations  -  one  for  the  interior  and  one  for  implementing  the  boundary 
conditions. 


4-1.  Interior  Formulation 

First,  consider  the  two-dimensional  wave  equation  (same  formulation  as  presented  in  (3)  with  U  =  V  =  0) 

-  clV2h  +  fh  =  0.  (7) 

Multiplying  by  the  test  functions  T,;  and  integrating  over  the  circular  domain  yields  the  weak  integral  form 


'lb 


d2h 

W 


dfl  — 


<F iX72h  <m  +  /  if2h  dn  =  0. 


Transferring  the  second  order  spatial  derivatives  from  h  to  the  basis  functions  via  integration  by  parts  and 
applying  the  divergence  theorem  to  recast  one  surface  integral  term  as  a  boundary  integral  gives  us 


, 0_  2 
dt2  dn  co 


Wifi  ■  V/t  dfl  +  cl  /  V'F;  •  Vh  dfl  +  /  ^ifh  dn  =  0. 


(8) 


Of  note  now  is  that  the  boundary  condition  (5)  contains  a  radial  derivatives  of  h  that  on  the  circle  is 
precisely  the  normal  derivative  n  ■  S7h.  This  allows  direct  implementation  of  the  boundary  condition  into 
(8)  as  follows: 


£  da-ell^  Vl/,:  if1  ~  "  h k )  ^  +  C°  L  ‘  V/l  ^  +  L  * if2h  ^  =  °'  (9) 

Here,  since  on  the  boundary  the  radius  is  fixed,  the  term  may  be  treated  as  a  constant. 

A  similar  weak  form  is  constructed  for  the  boundary  formulation  by  multiplying  (6)  by  the  test  functions 
Q  and  integrating  over  T  yielding  (after  by  integration  by  parts): 


1 

co 


dT+J- 
«  at  r 


Cifij 


dr 


dr 


1  c  1  end  1  f  d(i  d<t>j- 1 
4r2  1  d6  I  start  4 r2  Jr  80  d0 


dr. 


(10) 


We  now  use  the  fact  that  the  boundary  is  continuous  and  closed  to  surmise  that  the  endpoint  evaluation 
term  vanishes.  The  formal  problem  statement  is  then:  Find  h  G  V  and  <pj  €  Vr  where  j  =  1, . . .  J  —  1,  such 
that  Equations  (9)  and  (10)  are  satisfied  V  4b  £  V  and  Q  £  Vr- 


4 


4-2.  Galerkin  Expansion 

We  now  turn  our  attention  to  the  spatial  discretization.  First,  we  expand  the  solution  variables  h  and 
(j>j  using  the  same  basis  functions  used  in  the  weak  form  as  follows: 

Np  Nb 

hi v  =  ti >khk,  <i>jN = Tfc^ j ,  j  =  i,  2, . . . ,  j — l.  (it) 

fc=l  fc=l 

Here,  Np  refers  to  the  number  of  points  that  52  is  discretized  into  and  Nb  refers  to  the  number  of  points 
that  r £  is  discretized  into.  Next,  we  substitute  this  basis  function  expansion  directly  into  the  weak  form, 
resulting  in  the  following  matrix  form  of  this  problem: 


h  =  Ah,  +  Mh  +  C<f>i.  (12) 

If  we  examine  the  boundary  auxiliary  variable  formulation  ( 10)  we  see  that  the  selection  of  appropriate 
Cj  values  for  the  auxiliary  variables  has  not  yet  been  addressed. 

It  should  be  noted  that  van  Joolen  et  al.  [21]  show  that  any  choice  of  Cj  is  guaranteed  to  reduce  spurious 
reflection  as  the  order  of  the  NRBC  ( J)  increases.  While  we  omit  the  details  here,  the  core  of  this  argument 
is  the  computation  of  a  so-called  reflection  coefficient  that  is  a  product  of  J  factors,  each  of  which  are  less 
than  one.  The  reflection  caused  by  the  artificial  boundary  must  decrease  as  the  order  of  the  NRBC  increases. 
They  note,  “Of  course,  a  good  choice  for  the  Cj  would  lead  to  better  accuracy  with  a  lower  order  J,  but 
even  if  the  ‘wrong’  Cj ’s  are  taken  . . .  one  is  still  guaranteed  to  reduce  the  spurious  reflection  as  the  order  J 
increases.”  [21,  page  1045] 

If  we  now  collect  the  terms  on  the  left  and  right,  we  get  the  matrix  form  of  the  problem: 

E$  =  F$  +  h.  (13) 


4-3.  Time  Integration 

The  formulation  outlined  in  (12)  and  (13)  constitute  a  system  of  coupled  ODEs  that  must  be  solved 
to  yield  a  solution  for  h(x,y,t).  Since  the  goal  of  this  analysis  is  to  uncover  the  “true”  gains  made  by 
liigh-order  boundary  treatment,  it  is  possible  that  any  high-order  treatment  of  the  boundary  and  spatial 
discretization  without  considering  a  high-order  treatment  of  the  temporal  component  could  mask  gains  made 
by  a  high-order  boundary  treatment.  For  this  purpose,  our  approach  uses  standard  kth  order  Runge-Kutta 
(RK)  methods  (up  to  order  10)  to  integrate  the  system  in  time. 

The  set-up  of  this  scheme  is  a  standard  one,  namely,  the  second  order  system  is  expanded  to  a  larger 
system  of  first  order  ODEs,  then  solved  appropriately  using  the  associated  RK  tableau.  For  most  cases  in 
this  analysis  (unless  otherwise  stated)  time  integration  is  performed  using  a  4th  order  RK  scheme  using  a 
time-step  chosen  to  ensure  a  Courant  number  of  0.25,  where  the  Courant  number  is  defined: 

Courant  number  =  —  0 

y  (Ax)2  +  (Ay)2 

Here,  Ax  and  Ay  are  chosen  as  the  minimum  distance  between  any  two  points  in  the  x—  or  y—  directions 
respectively.  Additionally,  This  choice  is  made  since  the  interpolation  points  are  not  uniformly  distributed 
when  using  spectral  elements. 

4-4-  Results  for  the  HH  Formulation 

A  series  of  experiments  was  conducted  to  determine  the  effect  of  the  HH  boundary  condition  for  various 
SE  and  NRBC  orders.  Since  the  formulation  is  designed  for  circular  boundaries,  we  consider  only  circular 
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boundaries  with  unstructured  grids.  In  each  case,  we  choose  the  number  of  elements  to  yield  approxi¬ 
mately  3,  000  global  points.  We  take  a  smooth,  two-dimensional  “oval-shaped”  initial  condition  with  shape 
parameters  ax  =  g,  <Jy  =  3,  further  rotated  by  an  angle  of  9  =  The  initial  condition  used  here  is: 

h(x,y,  Q)  =  e-(a*2+2bxy+cy2),  h(x,y,  0)  =  0.  (14) 


Here,  the  parameters  a,  b ,  and  c  are  defined  as  follows: 


cos2  6  sin2  9  sin  2 9  sin  2 9  sin2  9  cos2  9 

+  2cry  4  o'2  +  4 0-2  c  2<t2  +  2cr2 


(15) 


Again,  the  solution  is  compared  to  one  computed  on  a  larger  domain  allowing  the  wave  to  propagate  out  of 
the  NRBC  domain  but  not  yet  impinge  on  the  non-physical  boundary  used  to  compute  the  solution  on  the 
larger  domain.  Qualitative  results  are  shown  in  Figure  4.7  and  quantitative  Lq  errors  are  shown  for  various 
NRBC  orders  for  SE  orders  up  to  6  in  Table  2.  No  further  improvement  was  observed  for  SE  orders  above 
order  6. 


j.5.  Adjustments  to  HH  to  Include  Mild  Dispersion 

The  unstructured  grid  (see  e.g.  Figure  3  from  [22])  representation  of  the  HH  formulation  has  been 
demonstrated  to  significantly  reduce  reflection  caused  by  the  boundary  for  the  standard  wave  equation.  The 
question  now  arises,  can  this  formulation  be  extended  to  include  dispersive  effects  such  as  Coriolis?  In  [17], 
van  Joolen  et  al.  presented  a  method  to  extend  the  HH  formulation  for  the  standard  wave  equation  under 
mild  dispersion.  While  this  formulation  was  well  grounded  mathematically  it  was  never  implemented.  A 
brief  synopsis  of  their  derivation  follows  with  results  presented  for  mild  dispersion  where  /2  =0.1. 

We  first  consider  the  KGE  without  advection  (in  polar  coordinates  as  in  the  HH  derivation): 


1  d2h_d2h  1  dh  1  d2h  f2 
Cq  dt2  dr 2  r  dr  +  r2  d92  Cg 


As  has  been  previously  discussed,  in  the  geophysical  context,  the  dispersion  parameter  is  typically  small. 
We  assume  here  that 


c0K 


<  1, 


(17) 


where  K  is  a  typical  wave  number  appearing  in  the  solution.  Now,  apply  the  Fourier  transform  to  (4)  and 
(16)  in  time  to  yield: 


d2h 
dr 2 
d2h 
dr2 


1  dh  1  d2h 
r  dr  +  r2  d92 
1  dh  1  d2h 
r  dr  +  r2  d92 


Wave 

Klein-Gordon 


where  w  is  the  frequency  and  h  is  the  frequency  domain  representation  of  h.  In  both  cases,  we  obtain  the 
Helmholtz  equation: 


K2h  + 


d2h  1  dh  1  d2h 
dr2  +  r  dr  +  r2  d92 


In  the  non-dispersive  case,  K  =  —  =  K  and  K  =  JK2  —  ^  in  the  dispersive  case.  In  order  to  facilitate 

Co  Y  Cq 

the  conversion  back  to  the  time  domain,  we  now  consider  a  Taylor  series  approximation  to  the  square  root 
term  found  in  the  dispersive  case,  i.e., 


\J  1  —  x  =  1  — 


\x  +  Q(x2). 
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Provided  that  x  is  small,  we  can  truncate  the  0(x2)  terms.  In  our  case,  from  (17)  we  can  reasonably  make 
this  assumption  yielding  for  the  dispersive  case: 


I  P  P  ^  „  P 

K  Y  CqK2  ~  1  2c2K2  2c2  A' ' 

We  now  see  that  in  the  frequency  domain,  an  equation  valid  in  the  non-dispersive  case  is  valid  in  the 
dispersive  case  if  we  make  the  replacement: 


/  f2  f2 

Jk2  -  4  «  K-  J 


2  c2K' 


(18) 


We  now  turn  our  attention  to  the  boundary  condition  (5)  and  (6)  that  we  Fourier  transform  in  time  to  yield: 

(19) 


l  dh  1  f  r 
—iKh  +  — — b  —  ft  —  <f>\ 
or  2  r 


Jl  {i~\)  X  1  32<A7-i 


4r2 


f'i-i 


4r2  dd2 


—  4>j+ii  j  1,  •  •  •  >  J  1- 


(20) 


Making  the  substitution  (18),  we  obtain  the  dispersive  version  of  the  HH  formulation  in  the  frequency 
domain,  i.e., 


—iKh  + 


?72 


o  ft  - b  ... —  4~  — ft  —  (b  i 
2cqAT  dr  2  r 


iP  7 

24K**  +  rl 


(j-h? 

4  r2 


Pj-i 


9ft.  1 

9r  2r 

1  d2(j)j- 1 
4r2  9ft2 


—  •Aj+i)  J  —  1)  ■  •  ■ ,  J  !• 


Transforming  these  equations  back  into  the  time  domain  results  in  the  final  HH  boundary  formulation 
for  the  KGE: 


pdh  ,  f 

cq  dt 


+  —  /  0j(r)  - ^2— 


+  bol  h^dT  +  ar  +  2rh  =  ,l‘1 

1  d2(t>j- 1 


m(t) 


Cq  9f  2co 


4r2  9ft2 


-  0j+i 


(21) 


(22) 


ft*) 


where 

j  =  1, . . . ,  J  —  1,  <p0  =2h  and  <f>j  =  0. 

It  should  be  noted  that  van  Joolen  et  al.  [17]  show  how  m(t)  and  n(t)  can  be  calculated  in  each  time-step  to 
keep  the  boundary  condition  local  in  time  without  having  to  store  and  operate  on  the  history  of  the  solution. 
For  this  analysis,  a  simple  trapezoidal  approximation  was  used  to  approximate  the  integral. 

The  weak  form  of  the  formulation  is  now  constructed.  We  consider  the  KGE  in  its  general  form: 


92ft 

~dP 


CgV2ft  +  ph  =  0. 


Multiplying  by  the  test  functions  ’F;  and  integrating  over  the  circular  domain  yields  the  weak  integral  form 

92ft 


1]ft  dn-c0 


dt,  V  ft  dfl  +  /  /  Tift  dH  =  0. 

:  J  o 
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Transferring  the  second  order  spatial  derivatives  from  h  to  the  basis  functions  via  integration  by  parts  and 
applying  the  divergence  theorem  to  recast  one  surface  integral  term  as  a  boundary  integral  gives  us 

[  [  %n-vh<Ki  +  <$  I  v^i-vhdn  +  f  [  Vihdn  =  o.  (23) 

Jo.  ot  J  p  Jq  Jq 

Of  note  now  is  that  the  boundary  condition  (21)  contains  a  radial  derivatives  of  h  that  on  the  circle  is 
precisely  the  normal  derivative  n  ■  S7h.  This  allows  direct  implementation  of  the  boundary  condition  into 
(23)  as  follows: 


t  nr,  2 

vIb  fjt2  dn~co 


+Co 


1  dh  1  f2  \ 

</>! - —  h-f—m(t))  dfl 

c0  dt  2 r  2c0  v  ’  J 

VT?;  •  Vh  dfl  +  f2  [  4>zh  dfl  =  0. 


(24) 


Here,  since  on  the  boundary  the  radius  is  fixed,  the  g-  term  may  be  treated  as  a  constant. 

A  similar  weak  form  is  constructed  for  the  boundary  formulation  by  multiplying  (22)  by  the  test  functions 
Q  and  integrating  over  T  yielding  (after  integration  by  parts): 


d(j)j 


f 


-  I  Ci-^  dT  +  /  Ci  n(t)  dT  +  j  dT  - 


co 


2  co 

_ 1  c  dfa-i 

4  r2  1  dd 


4  r2 


end  J 

start  4r2 


d( i  d(f>j- 1 


89  dd 


dT  = 


0  0j — i  dr 
r 

Ci4,j+ 1  dr. 


(25) 


We  now  use  the  fact  that  the  boundary  is  continuous  and  closed  to  surmise  that  the  endpoint  evaluation 
term  vanishes.  The  formal  problem  statement  is  then:  Find  h  £  V  and  < fij  £  Vr  where  j  =  1, . . .  J  —  1,  such 
that  Equations  (24)  and  (25)  are  satisfied  V  T,  £  V  and  Q  £  Vr- 


^.6.  Results  for  HH  with  Dispersion 

A  series  of  experiments  was  conducted  to  determine  the  effect  of  the  HH  boundary  condition  extended  to 
include  mild  dispersion  for  various  SE  and  NRBC  orders.  The  set-up  is  identical  to  the  experiments  without 
dispersion,  except  the  dispersion  parameter  is  set  to  / 2  =  0.1.  Qualitative  results  are  shown  in  Figure  4.7 
and  quantitative  L ^  errors  are  shown  for  various  NRBC  orders  for  SE  orders  up  to  6  in  Table  1.  As  in  the 
non-dispersive  case,  no  improvement  was  observed  for  SE  orders  above  order  6. 


4-7.  Effects  of  Time  Integration  Technique 

At  the  outset  of  this  work,  it  was  believed  that  at  some  point  the  improvements  realized  by  improving 
the  spatial  discretization  and  NRBC  would  eventually  be  limited  by  the  time  integration  scheme  [23].  To 
this  end,  the  order  of  the  time  integration  scheme  (RK2  -RK10)  was  varied  to  examine  the  effects  of  time 
integration  on  the  accuracy  of  the  solution.  As  has  already  been  presented,  gains  made  by  increasing  the 
order  of  the  NRBC  halt  for  lower  order  (second  order)  spectral  elements  after  J  =  5.  Even  for  high  order 
(order  8  and  16)  spectral  elements,  the  gains  made  by  increasing  the  order  of  the  NRBC  are  limited  at  some 
point  using  RK4.  In  [16]  the  authors  showed  that  high-order  time  integration  allowed  boundary  gains  to 
improve  solution  quality  for  the  KGE  under  zero  advection. 

For  this  experiment,  consider  the  KGE  on  a  semi-infinite  domain  with  h  =  0  on  T\y.  To  ensure  that 
any  boundary  or  time  effects  are  not  masked  by  the  interior  discretization,  we  consider  8th  order  spectral 
elements  discretized  into  9,409  global  points.  The  Gaussian  initial  condition  is  used  and  is  evaluated  until 
t  =  3.  The  reference  solution  in  this  case  was  computed  as  described  previously,  except  that  time  integration 
was  performed  with  a  10th  order  Runge-Kutta  scheme  (RK10)  using  a  time-step  chosen  to  ensure  a  Courant 
number  of  0.1.  Further,  to  exaggerate  the  issue  of  advection  and  dispersion,  we  consider  the  case  where 
U  =  0.75  and  /2  =  1. 


As  can  be  observed  in  Figure  5,  gains  made  by  improving  the  time  integration  matter  only  if  combined 
with  high-order  treatment  of  the  boundary.  Conversely  -  gains  using  high-order  treatment  of  the  boundary 
can  only  be  realized  if  there  is  a  sufficiently  high-order  treatment  of  the  time  integration.  While  this  analysis 
conducted  time  integration  using  up  to  RK10  (18  stages),  RK4  or  RK5  appears  to  be  sufficient  in  exploiting 
gains  when  using  higlr-order  (up  to  J  =  15)  boundary  treatment.  It  should  be  noted  that  these  results  (error 
on  the  order  of  10-9)  cannot  be  observed  unless  high-order  treatment  of  the  interior  also  accompanies  the 
high-order  treatment  of  the  boundary  and  time.  Several  experiments  were  conducted  which  varied  the  order 
of  the  interior,  boundary  and  time  integration.  The  clear  result  was  that  without  high-order  treatment  of 
all  components  in  concert,  convergence  to  the  reference  solution  is  stalled. 

While  these  results  are  for  a  specific  problem  (KGE  with  advection),  we  believe  that  the  principle  of  a 
balanced  approach  to  all  components  (interior,  boundary  and  time)  is  a  sound,  extensible  procedure  for 
any  problem.  If  high-order  treatment  of  any  of  the  three  components  is  missing,  the  high-order  treatment 
of  the  other  components  is  essentially  wasted;  the  results  reported  here  confirm  this. 


Acknowledgements 

The  first  author  is  indebted  to  the  US  Army  for  its  support.  The  authors  would  like  to  express  their 
appreciation  to  the  Naval  Postgraduate  School  for  its  support  of  this  research.  Finally,  the  authors  wish  to 
thank  the  reviewers  for  their  helpful  insights,  suggestions  and  comments. 

[1]  D.  Givoli,  Non-reflecting  boundary  conditions,  Journal  of  Computational  Physics  94  (1)  (1991)  1-29.  doi:10. 1016/0021- 
9991(91)90135-8. 

[2]  B.  Engquist,  A.  Majda,  Radiation  boundary  conditions  for  acoustic  and  elastic  calculations,  Communications  on  Pure  and 
Applied  Mathematics  32  (3)  (1979)  313-357.  doi:10.1002/cpa.3160320303. 

[3]  A.  Bayliss,  E.  Turkel,  Radiation  boundary  conditions  for  wave-like  equations,  Communications  on  Pure  and  Applied 
Mathematics  33  (6)  (1980)  707-725.  doi:10.1002/cpa.3160330603. 

[4]  F.  Collino,  High  order  absorbing  boundary  conditions  for  wave  propagation  models.  Straight  line  boundary  and  corner 
cases,  in:  R.  Kleinman,  et  al.  (Eds.),  Proceedings  of  the  2nd  International  Conference  on  Mathematical  &  Numerical 
Aspects  of  Wave  Propagation,  SIAM,  Delaware,  1993,  pp.  161  171. 

[5]  S.  V.  Tsynkov,  Numerical  solution  of  problems  on  unbounded  domains.  A  review,  Applied  Numerical  Mathematics  27  (4) 
(1998)  465-532.  doi:10.1016/S0168-9274(98)00025-7. 

[6]  T.  Hagstrom,  Radiation  boundary  conditions  for  the  numerical  simulation  of  waves,  Acta  Numerica  8  (-1)  (1999)  47-106. 
doi:10.1017/S0962492900002890. 

[7]  R.  L.  Higdon,  Absorbing  boundary  conditions  for  difference  approximations  to  the  multi-dimensional  wave  equation, 
Mathematics  of  Computation  47  (176)  (1986)  437—459. 

[8]  R.  L.  Higdon,  Radiation  boundary  conditions  for  dispersive  waves,  SIAM  Journal  on  Numerical  Analysis  31  (1)  (1994) 
64-100.  doi:10. 1137/0731004. 

[9]  D.  Givoli,  B.  Neta,  High-order  non-reflecting  boundary  conditions  for  dispersive  waves,  Wave  Motion  37  (3)  (2003)  257-271. 
doi:10.1016/S0165-2125(02)00074-4. 

[10]  D.  Givoli,  B.  Neta,  High-order  non-reflecting  boundary  scheme  for  time-dependent  waves,  Journal  of  Computational 
Physics  186  (1)  (2003)  24-46.  doi:10.1016/S0021-9991(03)00005-6. 

[11]  B.  Neta,  V.  van  Joolen,  J.  R.  Dea,  D.  Givoli,  Application  of  high-order  higdon  non-reflecting  boundary  conditions 
to  linear  shallow  water  models,  Communications  in  Numerical  Methods  in  Engineering  24  (11)  (2008)  1459—1466. 
doi:10.1002/cnm.l044. 

[12]  V.  J.  van  Joolen,  B.  Neta,  D.  Givoli,  High-order  boundary  conditions  for  linearized  shallow  water  equations  with 
stratification,  dispersion  and  advection,  International  Journal  for  Numerical  Methods  in  Fluids  46  (4)  (2004)  361—381. 
doi:10.1002/fld.754. 

[13]  J.  M.  Lindquist,  F.  X.  Giraldo,  B.  Neta,  Klein-Gordon  equation  with  advection  on  unbounded  domains  using  spectral 
elements  and  high-order  non-reflecting  boundary  conditions,  Applied  Mathematics  and  Computation  217  (2010)  2710- 
2723. 

[14]  D.  Givoli,  B.  Neta,  I.  Patlashenko,  Finite  element  analysis  of  time-dependent  semi-infinite  wave-guides  with  high- 
order  boundary  treatment,  International  Journal  for  Numerical  Methods  in  Engineering  58  (13)  (2003)  1955—1983. 
doi:10.1002/nme.842. 

[15]  L.  Kucherov,  D.  Givoli,  High-order  absorbing  boundary  conditions  incorporated  in  a  spectral  element  formulation,  Inter¬ 
national  Journal  for  Numerical  Methods  in  Biomedical  Engineeringdoi:  10.1002/cnm.ll88. 

[16]  J.  M.  Lindquist,  B.  Neta,  F.  X.  Giraldo,  A  spectral  element  solution  of  the  Klein-Gordon  equation  with  high-order 
treatment  of  time  and  non-reflecting  boundary,  Wave  Motion  47  (5)  (2010)  289-298.  doi:10.1016/j.wavemoti.2009.11.007. 

[17]  V.  van  Joolen,  D.  Givoli,  B.  Neta,  High-order  non-reflecting  boundary  conditions  for  dispersive  waves  in  cartesian,  cylin¬ 
drical  and  spherical  coordinate  systems,  International  Journal  of  Computational  Fluid  Dynamics  17  (4)  (2003)  263—274. 

[18]  V.  J.  van  Joolen,  Application  of  higdon  non-reflecting  boundary  conditions  to  shallow  water  models,  PhD  dissertation, 
Naval  Postgraduate  School,  Monterey,  California  (2003). 


9 


[19]  J.  Pedlosky,  Geophysical  Fluid  Dynamics,  2nd  Edition,  Springer- Verlag,  New  York,  1986. 

[20]  L.  L.  Thompson,  R.  Huan,  Accurate  radiation  boundary  conditions  for  the  two-dimensional  wave  equation  on  unbounded 
domains,  Computer  Methods  in  Applied  Mechanics  and  Engineering  191  (1)  (2001)  311-351. 

[21]  V.  J.  van  Joolen,  B.  Neta,  D.  Givoli,  High-order  higdon-like  boundary  conditions  for  exterior  transient  wave  problems, 
International  Journal  for  Numerical  Methods  in  Engineering  63  (7)  (2005)  1041-1068.  doi:10.1002/nme.l322. 

[22]  G.  Zhao,  M.  Xinwa,  Automesh  2-d:  A  fully  automatic  adaptive  quadrilateral  mesh  generator, 

http://www.automesh2d.com/default.htm  [Accessed:  February  12,  2010]. 

[23]  A.  R.  Curtis,  High-order  explicit  Runge-Kutta  formulae,  their  uses,  and  limitations,  IMA  Journal  of  Applied  Mathematics 
16  (1975)  35-53.  doi:10.1093/imamat/16.1.35. 


10 


Figure(s) 
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Figure  1:  The  shallow  water  model  with  irregular  bottom  topography.  (Adapted  from  [19,  page  58]). 
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Figure(s) 


(a)  Ref.  Solution  (t  =  1) 


(b)  Ref.  Solution  ( t  =  2) 


(c)  Ref.  Solution  ( t  =  3) 


(d)  NRBC  J  =  4  (t  =  1) 


(e)  NRBC  J  =  4  (t  =  2) 


(f)  NRBC  J  =  4  (t  =  3) 


Figure  2:  Open  Domain,  4™  order  spectral  elements  ( J  =  4)  using  oblique  Gaussian  initial  condition  shown  for  t  =  1,  2,  3.  Top 
Plots:  Contour  plots  of  reference  solution  solved  on  extended  domain.  Superimposed  black  circle  indicates  NRBC  domain. 
Bottom  Plots:  Contour  plots  of  various  NRBC  boundary  configurations  using  J  =  4. 
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Figure(s) 


Figure  3:  Sample  mesh  generated  using  Automesh-2D.  Geometry  for  the  letters  graciously  provided  by  [22]. 
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Figure(s) 


(a)  Ref.  Solution  (£  =  1) 


(b)  Ref.  Solution  ( t  =  2) 


(c)  Ref.  Solution  ( t  =  3) 


(d)  NRBC  J  =  4  (t  =  1) 


(e)  NRBC  J  =  4  (t  =  2) 


(f)  NRBC  J  =  4  (t  =  3) 


Figure  4:  Open  Domain,  4t/l  order  spectral  elements  (J  =  4)  using  oblique  Gaussian  initial  condition  shown  for  t  =  1,2,3 
under  dispersion  / 2  =  0.1.  Top  Plots:  Contour  plots  of  reference  solution  solved  on  extended  domain.  Superimposed  black 
circle  indicates  NRBC  domain.  Bottom  Plots:  Contour  plots  of  various  NRBC  boundary  configurations  using  J  =  4. 
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Figure(s) 


Figure  5:  Semi-infinite  channel  Lq  error  versus  RK  time  integration  order  using  various  order  NRBCs.  Domain  is  discretized 
into  9,409  points  for  8th  order  spectral  elements  with  advection  velocities  U  =  0.75,  V  =  0.0  and  dispersion  parameter  f2  =  1. 
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Table(s) 


Table  1 :  Lq  Error  as  a  function  of  NRBC  Order  for  Hagstrom  Hariharan  NRBC  formulation  using  various  spectral  element 
orders  on  the  circular  NRBC  domain.  Oblique  Gaussian  initial  condition  is  used  with  dispersion  parameter  set  to  /2  =0.1. 


NRBC  Order 

Error 

Linear  Elements 

Lq  Error 

Order  2  Elements 

L\ 2  Error 
Order  4  Elements 

Error 

Order  6  Elements 

J  =  1 

0.07290 

0.03555 

0.03369 

0.03293 

J  =  2 

0.02684 

0.00371 

0.00283 

0.00248 

J  =  3 

0.01869 

0.00258 

0.00192 

0.00204 

J  =  4 

0.01759 

0.00243 

0.00186 

0.00157 

J  =  5 

0.01744 

0.00240 

0.00181 

0.00154 

J  =  10 

0.01742 

0.00240 

0.00180 

0.00153 

J  =  20 

0.01742 

0.00240 

0.00180 

0.00153 
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Table(s) 


Table  2:  Lq  Error  as  a  function  of  NRBC  Order  for  Hagstrom  Hariharan  NRBC  formulation  using  various  spectral  element 
orders  on  the  circular  NRBC  domain.  Oblique  Gaussian  initial  condition  is  used. 


NRBC  Order 

Error 

Linear  Elements 

Lq  Error 

Order  2  Elements 

L\ 2  Error 
Order  4  Elements 

Error 

Order  6  Elements 

J  =  1 

0.09310 

0.04772 

0.04555 

0.04485 

J  =  2 

0.03381 

0.00465 

0.00355 

0.00315 

J  =  3 

0.02355 

0.00324 

0.00243 

0.00259 

J  =  4 

0.02217 

0.00305 

0.00236 

0.00201 

J  =  5 

0.02198 

0.00302 

0.00230 

0.00196 

J  =  10 

0.02195 

0.00302 

0.00228 

0.00196 

J  =  20 

0.02195 

0.00302 

0.00228 

0.00196 
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